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Abstract 


The paper is devoted to the description of the novel approach to calculate the vertical distribution of the 
average number of photon scattering events Ñ in a plane-parallel turbid layer having arbitrary local optical 
scattering characteristics and optical thickness. The derived equation for the vertical distribution of the 
average number of photon scattering events is a generalization of the well-known Ambartsumian’s formula 
(valid for homogeneous turbid media only) to the more general case of media with the vertically varying 
single scattering albedo. Numerical calculations are performed using the discrete ordinates technique in 
combination with the adjoint formulation of the radiative transfer equation. The performance of the 
numerical techniques introduced is studied using well-known asymptotic analytical equations valid for 
semi-infinite weakly absorbing layers. We also compare results obtained with published calculations of N 
derived from a superposition technique. 
© 2005 Elsevier Ltd. All rights reserved. 
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1. Introduction 


The information on the vertical distribution of photon scattering numbers M(t) (t is the optical 
depth counted from the top of a scattering layer) is of a great importance for a number of 
applications including the estimation of the photon horizontal transport in clouds and satellite 
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cloud remote sensing [1,2]. Historically, N(t) is calculated using the Monte-Carlo techniques that 
require a large computational time. Alternatively, a much more computationally efficient 
superposition technique was recently introduced [3]. 

This work is devoted to the presentation of a novel approach for the calculation of the vertical 
distribution of Ñ (t) using the discrete ordinates method (DOM) in combination with the adjoint 
formulation of the radiative transfer problem. The technique allows to obtain the values of N(z) 
with almost the same speed for layers of arbitrary optical thickness. 

The adjoint radiative transfer equation (ARTE) and its application to a number of 
environmental problems has recently been reviewed by Marchuk [4] and Box [5]. The idea 
behind this approach is the simultaneous solution of the standard (forward) radiative transfer 
equation (RTE) and also the adjoint RTE. Both solutions, as we show, give us the possibility to 
obtain the vertical profiles N(t). These profiles are related to the single scattering albedo profile as 
a function of optical depth and can be used in the information content studies in a number of 
inverse problems. 

The paper is organized as follows. In Section 2, we generalize the Ambartsumian’s formula for 
a case of vertically inhomogeneous turbid media. The Section 3 is devoted to the derivation of a 
simple analytical equation, which can be used to deduce the profile M(t) from precalculated 
diffused light intensities for the direct and adjoint RTEs. The Section 4 considers the application 
of this equation to a number of applied problems (e.g., for the calculation of N(t) for finite cloud 
layers sensed at multiple wavelengths). Here, we also show the comparison of our results with the 
exact superposition technique of Platnick [3] and the approximate analytical solution [6,7] valid 
for weakly absorbing semi-infinite media. A good agreement between these very different 
approaches is found. This indirectly confirms our calculations. 

Mathematical derivations related to the adjoint RTE are given in Appendix A. The discrete 
ordinates method (DOM) used in the numerical procedure is briefly discussed in Appendix B. In 
particular, we present the general solution of the homogeneous equation and the particular 
solution of the inhomogeneous direct and adjoint equations. 


2. The generalization of the Ambartsumian’s formula for the average number of photon scattering 
events 


2.1. Homogeneous layers 
The average number of photon scattering events in a random medium at the location specified 


by radius vector r and in the direction (0,ġ) can be found according to Ambartsumian [8] 
as 


N=——, (1) 
where œw is the single scattering albedo and / is the diffused light intensity at the considering 


location in the direction (0, @). Here 0 is the zenith angle and ¢ is the azimuth. The dependence of 
Ñ and J on variables r, 0, and ¢ is omitted for the sake of simplicity. 
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The proof of this formula is quite simple. Indeed, according to the original paper of 
Ambartsumian one can represent the intensity J as 


r=) ho, 2) 
n=0 


where Z, is the light intensity after n scattering events. 
One can find using Eq. (2) that the fraction P, of photons, which scattered exactly n-times is 
given by 
wo"L, 


sa a (3) 


This immediately gives Eq. (1) for Ñ = X% 9nPp. 

It is worth to notice that the same representation (2) for the scattered light intensity in turbid 
media without gaseous absorption has been used by van de Hulst [9] to calculate the mean photon 
path length /. Note that it follows for semi-infinite media: / = IÑ, where / is the photon mean free 
path. Also one can obtain for the average time of travel: 7 = //c. Here c is the speed of light in the 
medium. 

The main restriction of Eq. (1) is that œ has to be constant in the turbid medium under study. 
To generalize this formula for the case of w dependent on the vertical coordinate, we need to find 
yet another representation for N. For this purpose we represent the intensity as a series with 
respect to the variation of the single scattering albedo. Namely, we have in the vicinity of any 
point @: 


Iœ) = (6) + = 


(0-6), (4) 
@ 
where 4 | is the derivative of the function J with respect to œw calculated at the point @ and we 
neglect terms of higher order with respect to œ — æ. They should be accounted for if one needs to 
consider not only N but also variance, etc. 

Defining the variation of the intensity as AJ = (œ) — I(©) and the variation of the single 
scattering albedo as Aw = œ — ð, Eq. (4) can be rewritten in the form: 


Al(@) = a Aq. (5) 


w 


Dividing both sides of Eq. (5) by Z and dividing and multiplying right-hand side of this equation 
by œw we find 


L (6) 


(09) 


According to Eq. (6) Ñ can be treated as the proportionality coefficient between the relative 
variation of the single scattering albedo and the relative variation of the intensity at the point 
œ = @ in the framework of the linear approximation. 

Eq. (6) can be used to generalize Ambartsumian’s Eq. (1) for the case of œ dependent on t. Here 
t is the optical depth relative to the top of a light scattering layer. This problem is considered in 
the next subsection. 
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2.2. Inhomogeneous layers 


Let us assume that the single scattering albedo is a function of the optical depth, qt, i.e., w(t). 
Then a linear relationship between the relative variation of the parameters œ and J can be 
obtained by means of the Taylor series expansion of the intensity as a functional of the single 
scattering albedo. Namely, it follows: 


© 6I(a(t)) 
w(t) 


6a(t) dt, (7) 


@(t) 


Kot) = Ma(a)) + I 


where 6@(t) = w(t) — @(t), To is the optical thickness of a turbid layer, and 


To) _ imn TEO + 9a) — 1@@) 
da(t) — Ar>0 Jan w(t’) dT’ 


(8) 


is the variational derivative of J with respect to w(t). 

The variational derivative can be calculated, for example, with any radiative transfer model 
employing the numerical perturbation technique. In this case, the following approximation is used 
instead of exact Eq. (8): 


d1(@(t)) _ (a(t) + Aw@(tx)) — (a(t) — Aw(tx)) 
olt) 2A@(TK) 


(9) 


where Aw(t;) is the variation of the single scattering albedo profile at the level having optical 
depth t, and the radiation fields [(wm(t) + Aw(t,x)) are solutions of the RTE for two perturbed 
single scattering albedo profiles, w(t) + Aw(t,). Eq. (9) provides the approximation for the 
variational derivative of the intensity with respect to the parameter w(t) at the level t. By 
application of Eq. (9) at each discrete level, the derivative is constructed. 

Using the same transformations as those applied for the derivation of Eq. (6), we find that Eq. 
(7) can be rewritten in the form 


61(a(t)) 2 [ No) w(t) 

Kom) Jo c(t) 
where 6/(w(t)) = I(w(t)) — I(@(t)) and N(t) can be treated as the average number of photon 
scattering events in the infinitesimal layer positioned at the optical depth t. 

Summing up, we see that the formula for the calculation of the differential average number of 
photon scatterings is given by 
w(t) d1(o(t)) 

I 6a(t) ` 


dt, (10) 


N(t) = 


(11) 


If Ñ (t) is known for an arbitrary t, then the average number of photon scattering events in an 
inhomogeneous turbid medium in the layer positioned between t = 0 and tọ can be found as 


N= 4 N(a) dt. (12) 


Eq. (11) will be referred from now on as the generalized Ambartsumian’s formula. 
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The main problem for the practical usage of Eq. (11) is the calculation of the variational 
derivative. As it has been pointed out above, the numerical perturbation technique requires 
multiple solutions of the RTE. Therefore, our goal is to find an alternative expression for N(t), 
which does not include the variational derivative. In the next section we find the linear 
relationship between the variation of the intensity and the single scattering albedo without using 
Taylor expansions similar to Eq. (7). Instead we will consider the linearized form of the RTE. 


3. The linearized RTE and the alternative form of the generalized Ambartsumian’s formula 


We start from the standard RTE for a plane parallel light scattering layer illuminated by the 
radiation with the zenith angle 09 = arccos(u9). We assume further that: (i) the phase function is 
expanded in Legendre polynomial series, (ii) the light intensity is divided into the direct and 
diffuse parts, (iii) the diffuse light intensity is expanded in Fourier series, (iv) the underlying 
surface is the Lambertian one with the spherical albedo A. Under these assumptions we have the 
following form of the RTE for each Fourier harmonic m = 0,1,..., M (see [10] for details): 


dI” (t, u) 


(ea i / / t (Q) mn 
Pema was | Penbe a +5 oa (13) 
T 


with the boundary conditions 


T”(0,4)=0, u>0, 


1 
1°(t0, u) = Auge!" +2A f I(t, udp, =<, 
0 
I'"(t, u) =9, pw<0, m>0, (14) 


where t € [0, to] is the optical depth, to is the optical thickness of the medium, œ € [0, 1] is the 
single scattering albedo, u € [—1, 1] is the cosine of the polar angle as measured from the positive 
t-axis. It means that negative values of correspond to light propagated upwards; the opposite is 
true for downward fluxes. 

The scattering function p’"(c, u, u’) in Eq. (13) for the mth Fourier harmonic is 


K 
P(t, He!) = D> By, PRP), (15) 


k=m 


UL 


and the source function for the direct radiation Q” (t, u) is given by 


K 
O”, u) = D> Be PRWPR(ude (16) 


k=m 


where fy are expansion coefficients of the phase function in terms of Legendre polynomials. P¥ (u) 
are associated Legendre functions [11]. For the sake of simplicity we will suppress the explicit 
notation of the harmonic number m in further discussion. 
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It is convenient to rewrite Eq. (13) in the operator form 
w 
LI = 5 Olt.) (17) 
where Le is the linear differentialintegral operator, which combines all operations with intensity 
I(t, n) in Eq. (13) 
d œ f! p 
bean g+- | dpan) (18) 
dt 2 =] 


Taking into account that we need to derive an alternative expression for the variational 
derivative, we variate both sides of Eq. (17) and also the boundary conditions (see Eq. (14)) with 
respect to the single scattering albedo. The resulting operator equation for the variation of the 
intensity is written in the form: 


Lòl =% $.(2,1), (19) 


with the boundary conditions 


ôI(0,u)=0, p>0, 


1 
Sao = A f A E 20 20) 
0 


and the source function 
S(t, U) = Je(t, U) + Q(t, u) (21) 


with multiple and single scattering source functions defined as 
wo l £ 1 / 
J A(t, u) T 2 f vemen ` 
w 
Q) = 5 Q, (22) 


respectively. 

Eq. (19) provides the linear relationship between the variations of the intensity field and relative 
variations of the single scattering albedo under the assumption that the extinction coefficient 
remains constant. Thereby, this equation can be treated as the linearized RTE (LRTE) with 
respect to the variations of the single scattering albedo. 

The formal solution of the LRTE is easily found. Indeed, for the variation of the reflected 
radiance in the direction —y* at the level having the optical depth t* we have 


TO 
ôE = =| S(t, —u*)O(t = Tee 707 den(t) dt + (Wy, oD) 
Le Jo w(t) 


1 
+ 2Ae oT) * | I(t, we du. (23) 
0 
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Also it follows for the transmitted radiance in the direction p*: 


To N 
ôl = A f Se(t, UŽJOG* — rje DĂ eol dr + (Wt, ôT). (24) 
L® Jo co(T) 
Here ôI" = 61(t*, —u*), 6I' = ôI(1*, u*), and O(x) is the Heavyside step-function (see Appendix 
A) with the argument x = (t — t*). The usage of the Heavyside function allows us to extent 
integration limits over all range of t. 
For the sake of simplicity, we define the scalar product as 


(W4, 61) = f wan dt = i [ W(t, w)OT(t, K) du’ de, (25) 
where index g equals r for reflected and t for transmitted radiances, respectively, and 

Wile = FS pee, =, DO — ee (26) 

Wie = SS pe Wore — Deel 7) 


The first terms in Eqs. (23), (24) are the linear ones with respect to the variation of the single 
scattering albedo. But the second terms include the unknown variation of the intensity field. To 
exclude 6/ from this equation, we rewrite at first the last term in Eq. (23) in the form of the scalar 
product. Defining the function W;,(,7) as 

b(t, u) = 2Ae PVH" F(z — 9) UO(u), (28) 


where 6(t — to) and O(u) are the Dirac 6-function over t and the Heavyside step-function over u, 
respectively. 
Eqs. (23) and (24) can be combined in the following single equation: 


SII = i S4(z) d0(t) dt+(W!,61), g=r,t, (29) 
0 e c(t) 


where W(t, u) = Wi(t,u) + Wilt), W(t, = Wet, u), and 


1 ak x 
SE) == Ol — He" S(t, =u"), 
u 


1 4 ak x 
SO = za OG — De See, p) (30) 


To evaluate the scalar product containing ôZ on the right-hand side of Eq. (29), we use the well- 
known properties of the adjoint operator. By definition [12], for a given linear operator L, the 
adjoint operator L* satisfies the identity 


(g, Lf) = (L"g.f), (31) 


where f and g are two arbitrary functions from the domain of L. We assume that the function 
T (T, u) is the solution of the following linear operator equation: 


E = W”. (32) 
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Multiplying both sides of Eq. (32) by 6/ and using the definition of the adjoint operator (Eq. 
(31)), we have 


(ôI, W9) = (ôI, L*I*) = (LOI, I*) = (òS, I*), (33) 


where ôS is the right-hand side of the operator equation LôI = 6S. L contains all operations on ôT 
in the forward problem of radiative transfer (see [13] for details) and 6S = 6w/qS,(t, p) in our 
case. 

As it can be seen from Eq. (33), to find the alternative expression for the scalar product 
(W%, 651), we need to formulate the adjoint operator L* and find the solution 7 of Eq. (32) with 
appropriate boundary conditions. This problem is addressed in the Appendix A. The answer is 


w(t) 
w(t) 


where (,) is the integral over the u-variable as defined in Eq. (25). 
Therefore, it follows: 


(W4, 61) = [ (cai) dr, (34) 


a(t) 4 
olt i 
Comparing Eq. (35) with the expression (10) for the differential average number scatterings of 


photons, we find the generalized Ambartsumian’s formula (GAF), which does not include 
variational derivatives. It has the following form: 


ar = f ISOH (S15) (35) 
0 


Wa) = [SEC + (Selt), TO]. (36) 


[9 
The derived final formula is valid for the azimuthally averaged downward and upward 
propagated light fields at the optical depth t* in the direction +y*, where the positive sign is 
applied for the downward light filled. If the azimuthal dependence is required, one must find 
I(t, u) and I*(t, u) for all m = 0,1,..., M. Then the downward and upward light intensity are 
calculated as: 


1 M 
I(c*, tu", 6*) = a NOO = Som 0*, Eu") cosmo”), (37) 
m=0 


where ¢* is relative azimuth and 60 is the Kronecker delta. Finally, the expression for 
N(t,£u*, ġ*) is written in the form 


M 


Ñ’) = i NOSO) + (S20), FE,,(2))] cos(mg*), (38) 


m=0 


where the dependence of N%(t) and J’ on the u* and ¢* is omitted for the sake of simplicity. 

Eq. (38) is the main result of this work. It is worth pointing out that the derived formula is 
independent of the particular method used to solve direct and adjoint RTEs. The DOM is used for 
the solution of both direct and adjoint RTEs in this work. The brief discussion is given in 
Appendix B. 
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4. Results of numerical calculations 


Using Eqs. (82) and (83) (see Appendix B) for direct and adjoint radiances together with Eqs. 
(21), (22), (30) we have now all needed variables to calculate average number of photon scattering 
events according to the GAF (Eq. (36)) derived in the Section 3. We will consider here the 
azimuthally averaged case only and start with the case of a semi-infinite turbid layer. Then the 
solution can be found analytically under some assumptions (e.g., small light absorption in a turbid 
layer). 


4.1. Semi-infinite turbid layers 


The average number of scattering events for photons reflected from a semi-infinite weakly 


absorbing turbid layer is given by the following approximate equation [6,7] 
_ 2u 
=T. (39) 


where k is the diffusion constant of the radiative transfer theory [9] and 
_ K(u)K(u) 
u = ———_, 
R(Uo, 4, Q) 


The function R(t, 4, @) is the reflection function for the semi-infinite nonabsorbing layer [14,15]. 
It determines K(uọ) via the following relationship ([14)]): 


(40) 


1 
K(uo) = 0.75 (u +2 f MA TAD) dn), (41) 


where r(uo, 4) is the azimuthally averaged reflection function of a semi-infinite nonabsorbing layer 


1 27 
Hii =. [Rodded (42) 


Let us check the results of calculations using Eq. (36) against Eq. (39). This will prove Eq. (36) 
in the most difficult for calculations case involving large values of t and œ close to 1. Clearly, we 
have 


n= f Fod (43) 


The results of calculations of A for different values of the single scattering albedo as functions of 
the solar angle are given in Fig. 1. The well-known Heney—Greenstein (HG) phase function has 
been used in these calculations with the asymmetry parameter g equal to 0.85. This value of g is 
common for water clouds. Both results obtained from Eq. (36) and the approximate Eq. (39) are 
shown. Functions R(uọ, 4, ġ), K(u) and the constant k (see Eq. (39)) are taken from tables 
prepared by Yanovitskij [16] using exact radiative transfer calculations. We see that Eqs. (36), (43) 
produce the same results as the asymptotic formula (see Eq. (39)) at œ = 0.9999. The relative 
difference between these different approaches to calculate ñ is less than ~ 0.2% (see Fig. 1, right 
panel) in the broad range of solar zenith angles. The differences are larger for the smaller values of 
œ (larger absorption). As it can be deduced from Fig. 1, the maximal difference between two 
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Fig. 1. The integrated average numbers of scatterings for reflected photons for the case of a semi-infinite homogeneous 
cloud as a function of the cosine of the solar zenith angle (lines—numeric results, symbols—asymptotic results) at the 
optical thickness 1000, the Heney—Greenstein phase function with the asymmetry parameter equal to 0.85, single 
scattering albedos equal to 0.999 and 0.9999, and the cosine of viewing zenith angle u = 1. The right panel gives the 
relative error 6 = 1 —n,/n, of Eq. (39) in percent. Here na corresponds to the asymptotical result and ne is obtained by 
the numerical technique introduced in this paper. 


formulae reaches ~ —1.5% for œw = 0.999 and shows the pronounced dependence on the solar 
zenith angle. It confirms that the asymptotic formula (39) can be used for weakly absorbing media 
only. It is interesting to note that the accuracy of Eq. (39) is the best at the solar angle 60 degrees 
(see Fig. 1, right panel). For this particular geometry, Eq. (39) is valid for arbitrary w, which is a 
peculiar result taking into account the assumptions used in the derivation of the asymptotic 
solution [7]. Numerical calculations support that ñ ~ (1 — w) 2. Such a dependence was also 
found by Uesugi and Irvine [17]. Also we see that ñ increases with u. This could be easily 
understood on physical grounds. Indeed, the photons injected along the vertical in the 
nonabsorbing semi-infinite turbid medium with forward-extended phase functions need more 
time (and larger number of scatterings) to escape from the medium in the nadir direction than 
photons incident at grazing angles (small values of u, see Fig. 1). 

It is of interest to establish the range of the validity of Eq. (39). We adopt further the following 
approximate equations for R(uọ, 1,0) and K(uo) [18]: 


3 
K(Ho) = 7 (l + 2m), (44) 


0.37 + 1.94uọ p(n — arccos(uo)) 
1+ uo 4(1 + uo) Í 


R(uo, 1,0) = (45) 


where p(@) is the phase function. The constant k was assumed to be equal to J 3(1 — w)(1 — œg), 
which is a good approximation for values of w close to one [9]. Here g is the asymmetry 
parameter. The accuracy of these approximations has been investigated by Kokhanovsky [14,15]. 
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Fig. 2. The same as in Fig. 1 but as the function of the solar angle for the single scattering albedos equal to 0.9, 0.99, 
0.999, and 0.9999. 


The results of calculations of ñ for four different values of œ using described above 
approximations for related functions is shown in Fig. 2. The HG phase function has been used in 
calculations. At first we compare the results for œw = 0.9999 and œ = 0.999, which have been 
presented in Fig. 1 without approximations given by Eqs. (44), (45). We see that the relative 
difference between numerical and asymptotical values of ñ is larger now as it should be. So for 
œ = 0.9999 the maximal error is ~ 2% in the solar angles range 10—70°. It reaches ~ 6% in the 
vicinity of the solar angle ~ 80°. The error enhancement due to the usage of the approximations 
given by Eqs. (44), (45) for related functions is not so pronounced at œ = 0.999. In this case the 
maximal error is ~ 4% in the solar angles range [10—70°]. It reaches ~ 6% in the vicinity of the 
solar angle ~ 80°. 

The increase in absorption (œw = 0.99 and œw = 0.9) results in further increasing the relative error 
of Eq. (39) (especially for small solar angles and in vicinity of the solar angle ~ 40°). As can be 
seen in Fig. 2 the maximal error in this case can reach ~ 13% and ~ 40% for œw = 0.99 and 0.9, 
respectively. It is interesting that there is a minimum in the dependence of the average number of 
scatterings on the solar angle for media with appreciable absorption (see Fig. 2 at œ = 0.9). The 
minimum is located approximately at the solar zenith angle 09 = 40°. Such a minimum does not 
exist for values of the single scattering albedo very close to one (see Fig. 1). 

Concluding we can say that the asymptotic formula (Eq. (39)) together with approximations for 
related functions (see Eqs. (44), (45)) has the accuracy better than ~ 5% for nadir observation in 
broad range of the solar angles and w>0.999. For larger absorption, the accuracy is reduced 
considerably. However, the relative error is smaller than ~ 5% even for values of the single 
scattering albedo ~ 0.9 for solar angles 60—75° and the nadir observation. 


4.2. Comparison with the superposition technique 


Until now we have considered only results for the integrated average number of scatterings for 
reflected photons (ñ). In this section we compare results for the vertical distribution of the average 
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number of scatterings per the differential layer for reflected and transmitted photons (N'(t) and 
N t), respectively), calculated according to the generalized Ambartsumian’s formula (Eq. (36)) 
with results obtained using Platnick’s superposition technique (ST) [3]. 

The comparison of both results in four spectral bands usually used in cloud remote sensing 
problems is shown in Fig. 3. The results were obtained using the water cloud model having an 
optical thickness tọ = 8 and an effective cloud droplet radius 10 um. Mie theory was used to find 
local optical characteristics of a cloudy medium at wavelengths specified in Fig. 3 [2]. The cosines 
of the incident and observation angles were equal to 0.65 and 0.85, respectively. Only the 
azimuthally averaged values are shown in Fig. 3. We see that two different numerical techniques 
produce very similar results. This confirms their validity. Note that the ST has been compared 
with Monte-Carlo (MC) results as well, and excellent agreement was found (details in [3]). 

It follows from the analysis of Fig. 3 that functions N'(t) for reflected light are much more 
peaked. The positions of peaks correspond to layers from which the maximal contribution to the 
reflected light intensity comes. This depth generally decreases with the wavelength as shown in 
Fig. 3. Such peaks do not exist for transmitted photons. So contributions from different levels in 
the cloud are approximately equal for the transmitted light. It also means that droplet effective 
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Fig. 3. Vertical distribution of the average number of scatterings per differential layer for reflected (R) and transmitted 
(T) photons at wavelengths 0.66, 1.6, 2.2, and 3.7 um calculated for a homogeneous cloud with optical thickness 8 and 
the droplet effective radius 10 um. The cosine of the solar and viewing zenith angles were 0.65 and 0.85, respectively 
(symbols are calculations according to Platnick [3], lines are our technique). The results are averaged with respect to the 
azimuth. 
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radii in vertically inhomogeneous clouds sensed by passive reflection and transmission techniques 
will differ. 

Data similar to those shown in Fig. 3 are presented in Fig. 4 for A = 1.6 um but for the cloud 
optical thickness equal to 10 and for the nadir observation. The incidence angle was assumed to be 
equal 60°. Different numbers in Fig. 4 corresponds to the different microphysics of clouds. 
Namely cases of homogeneous clouds with effective radii 4, 10, and 16 um are shown by lines. The 
results for the vertically inhomogeneous layer are given by stars. To model the influence of the 
vertically inhomogeneity we assume here that the effective radius decreases linearly with 
the vertical coordinate z from 16 um at the top of a cloud to 441m at the bottom. We see that the 
dependence N'(t) for a vertically inhomogeneous cloud with values of the effective radius in the 
range 4-16 um can be quite well represented by the value of the average effective radius aef = 
10 um. 

It follows from Fig. 4 that the number of scattering events has a maximum at the optical 
thickness t* = 0.75—-1.0 counted from the cloud top position at the wavelength 1.6 um. 

In Fig. 5, similar results are shown for the wavelength 2.2m. In this case with larger 
absorption the value of t* is in the range 0.5—0.75. It means that different wavelengths are 
sensitive to different layers inside the cloud. 

In particular, the size of droplets usually increases with the height. This means that 
measurements at more absorbing wavelengths will give larger retrieved values of the 
cloud effective radius. This was noted by Platnick as well [1]. Sizes of ice crystals tend to 
decrease with height, so an opposite tendency may be expected for ice clouds away from 
nucleation regions. 
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Fig. 4. Vertical distribution of the average number of scatterings per differential layer for reflected photons calculated 
for homogeneous and inhomogeneous clouds with the optical thickness 10 and cosine of solar and viewing angles of 
Ho = 0.5 and u= 1, respectively. 1—effective radius 4, 2—10 um, 3—16 um, 4—vertically inhomogeneous cloud as 
described in the text. 
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Fig. 5. The same as in Fig. 4 but for 2 = 2.2 um. 


4.3. The vertical distribution of the average number of photon scattering events in multilayered cloud 
atmosphere with gaseous absorption 


In this section, we consider results for the vertical distribution of the average number of 
scatterings for a more realistic case of the atmosphere with Rayleigh scattering, aerosol scattering 
and absorption, gaseous absorption and multilayered clouds. Calculations have been carried out 
in the spectral range of the O2—A absorption band which is often used for the determination of the 
cloud top height (see for details Rozanov and Kokhanovsky [19] and references therein). We have 
used here the same model of atmosphere for pressure and temperature profiles, gaseous 
components, and aerosol optical parameters profiles as in our previous paper [20]. The cloud 
optical parameters were obtained for the Deirmendjian’s Cloud C.1 droplet distribution with 
effective radius 6um [14]. We discuss here only results for reflected photons at the top of 
atmosphere (60 km) in the nadir direction for the solar zenith angle 60°. This allows us to ignore 
azimuthal dependencies in the generalized Ambartsumian formula (Eq. (36)). 

In order to calculate the average number of scattering events in the whole atmosphere we 
present Eq. (12) in the more convenient form 


‘melt To n L e 
Ñ= f N(t)dt = 2 Ñ;, (46) 


where L is the number of layers in the whole atmosphere, and Ñ; is the differential number of 
photon scatterings per a finite layer 


TH 
Ñ; = N(t) dr, (47) 


qj 


with t; and t+; as optical depths at the top and the bottom of each layer, respectively. 
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The vertical distribution N; (left panel) for the atmosphere containing two clouds located 
between 2 and 4km (layer 1) and 7-9 km (layer 2), respectively, is shown in Fig. 6. In the right 
panel of this figure, the vertical profiles of the single scattering albedo used in the calculations are 
given. The integration in Eq. (47) has been done over 0.1 km sub-layers in clouds and over | km 
sub-layers in the free atmosphere. The optical thickness of the whole cloud system is constant for 
all cases shown in Fig. 6 and equals 20. 

The analysis of data given in this figure allows us to reach following conclusions: 


1. The vertical profile M(t) per a finite layer strongly depends on the distribution of the optical 
thickness between cloud layers in the presence of gaseous absorption. 

2. In each cloud layer the integrated value of the average number of photon scattering events 
strongly depends on the optical thickness and increases with the optical thickness. 

3. In the case of coinciding optical thicknesses of the upper and lower cloud layers, the average 
number of photon scattering events in the upper layer is always greater than that in the lower 
one. This can be explained taking into account that the single scattering albedo in the upper 
cloud layer is larger than that in the lower layer due to decreasing of the oxygen absorption 
coefficient for higher levels in a cloudy atmosphere. 

4. The sum of the average number of photon scatterings in the both cloud layers is smaller than 
average number of photon scatterings in the whole atmosphere with account for the molecular 
and aerosol scattering. This is due to a small contribution of the additional photon scattering 
on molecules and aerosol particles in the cloudless atmosphere. 
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Fig. 6. Vertical distribution of the average number of scatterings per finite layer for reflected photons at the top of 
atmosphere in the O.—A absorption band (A = 763.525 nm): Lines 1—4 are for a two-layered cloud system with optical 
thickness as follows: I—t, = 19, t = 1, (N’ = 23.4); 2—1 = 15, t2 =5,(N = 21.8); 31 = 10, t2 = 10, (WN = 23.4); 
4—1, = 5,1) = 15,(N = 26.2); 5—a single cloud (the top is located at 9km and the bottom is located at 2km) t = 20, 
(WN = 22.5). Numbers in brackets give the average number of scatterings for each case. The optical thickness of the 
oxygen layer is equal to 0.26 in the case under consideration. 
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The calculations of the vertical distributions Ñ; at the wavelengths of O2—A band with larger 
absorption show qualitatively the same dependencies as given above. However, the values of N; 
are smaller due to larger absorption (smaller values of the single scattering albedo). 


5. Conclusions 


We derived here the generalized Ambartsumian’s formula for the average number of photon 
scattering events in a turbid layer having arbitrary optical properties, which can vary along the 
vertical direction. The final equation is expressed via solutions of the direct and adjoint RTEs. 

Results obtained have been applied to the solution of selected raditaive transfer problems (see Figs. 
1—6). Our derivations and numerical calculations are of importance for atmospheric inverse problems 
as well. This is due to the fact that functions Ñ (t) coincide with weighting functions with respect to the 
vertically varying single scattering albedo. Also the value of Ñ allows to fund the average optical path 
length in a scattering medium with t —> oo. The root-mean-square horizontal photon optical path 
transport of reflected and transmitted photons can be easily found if N is known. 

The described approach for the calculation of the vertical distribution of the average number of 
photon scattering events is implemented in the software package SCIATRAN 2.0. The 
SCIATRAN 2.0 is freely available for non-commercial use at the website www.iup.physik.uni- 
bremen.de/sciatran. 
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Appendix A 
A.l. The adjoint radiative transfer operator 


To derive the adjoint operator in Eq. (32), we follow the approach suggested by Ustinov [13]. 
At first we rewrite the boundary conditions (Eq. (14)) of the direct RTE in the operator form 
Ld = St, 
Lgl = Sp, (48) 
where Są and Sp can be in the fact arbitrary but in our case (see Eq. (14)) St = 0 and S, = 
Apye—/“, Boundary conditions operators L; and Ly are defined as 


pa f * òO, 


To 1 
Lg = / dtd(t — Tt) low & —2A l. amua] ; (49) 
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where 0(t), 6(t — to) are Dirac 6-functions, A(u) = uO(u), O(u)—heavyside step-function: 


1, u>0, 
OW =) 9 neo (50) 


and symbol & is used to denote an integral operator, not a finite integral. Multiplying now the 
first equation in (48) by ô(t)A(u) and second one by 6(t — t)A(—p) and adding both these 
equations with Eq. (17) we have 
[Le + 0(t)A(W) Lt + O(t — To)A(—p) Lol 
= Qe + d(t)A(u)St + ÔC — T0)A(— u) Sp. (51) 
Eq. (51) is the generalized form of the direct RTE, which is equivalent to Eq. (17) but contains all 


operations with the radiance field including boundary conditions. According to Eq. (51) the 
generalized form of the direct operator can be written as 


L= Le + d(t)A(uyLy + ÔC — to)A(—p) Lo. (52) 


To find the adjoint operator we can now directly apply the definition of the adjoint operator 
according to Eq. (31). We can write after some algebra (see also [13]): 


L* = Le + ÒM LE + ÒC — 10) A(W)L5, (53) 
where L* is 


d 


1 
@ / / 
dt 2 —1 


The explicit expressions for L¥ and L% are 


= f "Oane: 
0 
1 


i= f ' dtd(t — To) low ® -24 | TETI : (55) 


Since the adjoint operator L* is defined, we can use now Eq. (32) to formulate the generalized 
form of the adjoint RTE which has to be solved to find I * needed for Eq. (33). 
Inserting function W(t, u) according to its definition after Eq. (29) we have: 


L*I* = W*(t, p) + (t — to) A(u)2AC P/E 
i = Wit.) (56) 


where Wi(z, u) and Wt(t, u) are defined according to Eqs. (26) and (27), respectively. 
Comparing right-hand side of Eqs. (56) with Eq. (53) we can rewrite both equations in the 
normal form of the RTE. Namely, we have for the case of reflected light: 


Le TF = Wit, u) (57) 
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with the boundary conditions 
r(0,u) =0, <0, 


0 
T¥(to, u) = 24e — 2.4 f I¥(to u)udu, p>0, (58) 


and for the case of transmitted light 
LÈT = WiC, p) (59) 
with the boundary conditions 
T*(0, w) = 0, u<0, 
0 
Fo, =—24 [onde u>0. (60) 


If solutions of Eqs. (57), (59) are found, then the scalar product (33) can be expressed in the next 


form 
Tena To z 0@(T) 
(W!,61) = f (Sul) E 


where the right-hand side of Eq. (19) have been used in Eq. (33) instead of 6S. 


dr, g= hr (61) 


Appendix B. The discrete-ordinates solution of the direct and adjoint RTEs 


To calculate M(t) according to Eq. (36) we need to find the solution of direct and adjoint RTEs. 
For this purpose we adopt here the DOM suggested by Chandrasekhar [21] (namely, its numerical 
implementation given by Stamnes et al. [22]). It allows us to find the general solution of the 
correspondent homogeneous equation. To find a particular solution of the inhomogeneous 
equation, the infinite-medium Green’s function approach suggested by Case and Zweifel [23] in 
the suitable for DOM form presented by Siewert [24] is used. 

According to the DOM, the radiance field is divided into N up-welling and N down-welling 
streams, producing radiance vector pairs I_(t) and I,(t), which are defined as 


I(t) = U(t, +u), I(t, +u), EE I(t, +u), (62) 


where u; are the quadrature points of the double-Gauss scheme adopted here. The detailed 
description of this scheme is given, for example, by Thomas and Stamnes [10]. The symbol T 
means the transpose vector. 

The general solution of the discrete-ordinates form of the RTE equation can be represented as 


L0) = $0) + EC), (63) 


where I (t) and IÈ (t) are the general solution of the homogeneous equation and a particular 
solution of the inhomogeneous equation, respectively. 

We assume further that inhomogeneous medium is divided in L homogeneous layers. The /th 
layer has the optical coordinate (depth) t;_; at the top and, respectively, t; at the bottom. The 
solution (see Eq. (63)) has to be found in all layers. 
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B.1. The general solution of the homogeneous equation 


The solution of the homogeneous equation in each homogeneous layer can be found in the form 


h(t) = Ove, (64) 


where v is the arbitrary constant and ®.(v) are two N-component vectors 


(v) = [dQ, +u), ply, +u), s oi, ey (65) 


The vector Dve Y can be treated as an elementary solution of the homogeneous RTE [25]. 
Taking into account that ®,(—v;) = ®_(v,) (see for details [24]) we see that the vector D_(v,Je*/ ‘i 
gives yet another elementary solution of this equation. The elementary solution remains a solution 
of the studied equation if we multiply it by an arbitrary constant. Summing up, we find the 
generalized solution of the homogeneous equation in the form 


N 
È (1) = `> [4t D (ve OY ais BO (v) E], (66) 
j=l 


where vectors ®(v;) and separation constants v; can be found after the solution of the 
appropriate eigenvalue problem [24]. 

For further simplification we introduce the 2N x 2N matrix ® as 
| ®+(V1),---, P (vw), O_(1),-.-, BC) 
~ | @_(1),...,®8_(vw), D40)... By (vy) |? 


D (67) 


the vector H containing 2N arbitrary constants 

H = [4}, A}, ey A B®, BB, , BA, (68) 
and 2N vector of the intensity field I'(z) = Gol. Under this assumption the short 
expression of the homogeneous solution Eq. (66) is 

I(t) = @E(1)H, (69) 


where E(t) is the 2N x 2N diagonal matrix, whose N first elements are e*-™-/" and N last 
elements are e7679, respectively, and j = 1,2,...,N. 


B.2. The particular solution of the inhomogeneous equation 


The infinite-medium Green’s function for a general form of the discrete-ordinates approxima- 
tion to the transport equation in plane geometry has been constructed by Barichello et al. [25]. 
The Green’s function is then used to define a particular solution of the inhomogeneous version of 
the discrete-ordinates equation. Taking into account the symmetry of the double-Gauss scheme 
used here, we will adopt a special version of this particular solution. According to Siewert [24] we 
write the particular solution in the adopted above vector—matrix form as 


P(t) = P(t), (70) 
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where the vector—function P(t) consists of 2N t-depending components. Namely, we have 


P(t) = [A] (©), A5(z), a , AÑ (T), BY (1), B5 (1), is Bol, (71) 


where components A?(t) and B?(z) are defined as 
AP) = C(t = ti : Vj, He Aa., 
BY(t) = S(t; — T : Vj, uo)e 77 b?, (72) 
p 


and functions C(t: x,y), S(t: x,y) and constants a} bP are defined according to Siewert [24]. 
Having defined a particular solution, we can now rewrite Eq. (63) in the form 


I(t) = ®E(t)H + ®P(7). (73) 


This formula gives the solution of the direct RTE for each homogeneous layer. The solution 
found includes the vector H composed of 2N arbitrary constants. H can be found using the 
boundary conditions and also the continuity of the intensity across layer interfaces. 


B.3. The solution of the adjoint RTE 


The direct and adjoint RTEs are closely related. Thereby, as has been pointed out in [26-29], 
the computational methods for the adjoint RTE can be developed by using the computational 
methods known for the direct RTE. Following [27], we state that the ARTE (Eq. (32)) can be 
written in the form of the direct RTE for the ‘primed’ intensity as 


LI (t, u) = W(t, —L), (74) 


where L, is the direct radiative transfer operator (see Eq. (18)) and the index g equals to t or r for 
the transmitted or reflected light case, respectively. Taking into account the explicit form of 
W(t, —p) given by Eqs. (26) and (27) we can conclude that the derived RTE for the ‘primed’ 
intensity is equivalent to the RTE with the so called upward and downward parallel surface source 
(PSS), which has been discussed by Qin et al. [30]. 

For the solution of Eq. (74) one can use the discrete ordinates approach. Then we find that Eq. 
(74) in the discrete-ordinates form has the same elementary solutions as those already described 
above. So the general solution for the ‘primed’ intensity can be written as 


I(t) = @E(2)H* + @P*(x), g= tr, (75) 


where Hy is the 2N-vector of the arbitrary constants. Components of the vector P” can be derived 
according to [24,25] after the integration of the right-hand side of Eq. (74) with the infinite- 
medium Green’s function. For the sake of simplicity we will assume here that t* (see Eq. (23)) 
coincides with the one of interfaces. Then we have for components of the vector P% after the 
integration: 

Reflected radiance (upward PSS): 


e = _ ak) / pk 
ABEC) = C(t — ta 2 vj WHET ape, 


BP“ (1) = S(t — t : vj, phen BR, (76) 
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Transmitted radiance (downward PSS): 
AB (2) = Slo = tra 2 ype“ OO ab, 


BE (2) = Clay — 0: vj ue 8 Be (77) 
where 


l 
= [O! (v)WW4 + @' (WW), g=r,t, 


*) NO) 
g 1 T 9 T g 
Bee = xO) [© (0) WWs + OT, WW], g=r,t (78) 


and 
WI = [W9 (+m), W9 (£m), ..., Wty), g=r,t, 


f Ol 
W"(+4) = ae Pi HŽ, =u), 


(69) 
Wey) = zp Palit EM). (79) 


Functions C(t : x,y), S(t: x,y), and N(v;) are the same as used for the particular solution of the 
direct RTE, W is diagonal matrix whose elements are quadrature weights (see for further details [24]). 

In Eq. (79) œ; and piu", +u) are single scattering albedos and phase functions, which are 
assumed to be constant inside each layer but can vary from one layer to another. 

Eq. (76) is used for layers with the optical depth larger than t*. Otherwise, we have: Ag C) = 


B%;(©) = 0. Eq. (77) is used for layers with the optical depth less than 1* (and ARM (t) = BR Ao) = 
0, otherwise). 

Having defined the solution for the ‘primed’ intensity, we can use now the inverse 
transformation /*(t, u) = T'(t,—u) to have the desired adjoint solution [5,26,27]. Employing such 
transformations to the vector solution of Eq. (75) we obtain the discrete-ordinates solution of the 
ARTE: 


1G) = TOEC)H} + TOP} (0). (80) 
Here 2N x 2N matrix T has the form 
T vet 81 


where 0 and 1 are N x N zero and diagonal identity matrices, respectively. 


B.4. The complete solution 


The complete direct and adjoint solution in the /th layer can be now written as 


V(t) = @E'(.)H! + ©'P'(0), (82) 
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I(t) = TO'E!(7)H*! + TO'P), (83) 


where / = 1,2,..., L and g =t,r. We need to find 2N x L components of the constant vectors H’ 
and HY . Using boundary conditions only 2N components are found. Remaining 2N x (L — 1) 
components have to be found using the continuity of the solution across layer interfaces. Detailed 
description of the appropriate matrix formulation is given by Thomas and Stamnes [10]. The 
explicit expression for the matrix of the system of linear equations with the vector—matrix 
notation in the form very close to one adopted here is given by Qin et al. [30]. 

It is worth to notice that the linear system matrix has the bound-structure and is independent of 
the sources. Thereby, for the numerical solution we have used here very efficient subroutines 
DGBTREF and DGBTRS from the LAPACK collection [31]. 
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